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A.  ACCOMPLISHMENTS 


The  following  is  a  brief  summary  of  accomplishments  under  Grant  AFOSR-84-0131.  Reports  and  papers 
listed  in  Section  B  contain  more  detailed  presentations  and  may  be  obtained  from  the  investigators. 


1.  The  Reduced  Basis  Method  for  Algebraic  Systems 

1.1  Subspace-Projector  Pairings.  Implementation  of  the  reduced  basis  method  requires  the  choice  of  a  subspace 
and  a  projector  onto  that  subspace.  For  an  arbitrarily  chosen  subspace-projector  pair,  existence  of  the  true  solu¬ 
tion  curve  is  not  sufficient  to  guarantee  the  existence  of  the  corresponding  reduced  basis  solution  curve.  How¬ 
ever,  when  the  former  curve  exists,  it  has  been  shown  in  [Al]  that  there  are  infinitely  many  subspace-projector 
pairings,  each  utilizing  an  arbitrarily  selected  subspace,  under  which  the  reduced  basis  solution  curve  exists. 
Moreover,  the  resulting  error  estimates  are  of  the  same  nature  as  those  that  apply  in  the  more  familiar  case  when 
a  subspace  is  paired  with  its  orthogonal  projector. 

12  Reduced  Bads  Additive  Correction  Methods.  Additive  Correction  Methods  have  been  considered  by  a 
number  of  authors  as  a  means  of  accelerating  slowly  converging  iterative  processes  (see,  for  example,  [A2]- 
[A4]).  Furthermore,  it  has  been  recognized  that  additive  correction  is  central  to  the  basic  idea  of  multigrid 
methods  [A4-A3].  Although  the  reduced  basis  method  in  its  original  form  appears  to  have  little  in  common 
with  additive  correction,  a  class  of  such  methods  has  been  developed  using  the  reduced  basis  concept  [1]. 
Furthermore,  it  has  been  shown  that  in  their  "two-grid"  form,  certain  multigrid  methods  are  special  cases  of  this 
class.  The  reduced  basis  point  of  view  provides  insight  into  the  error  reduction  capability  of  such  multigrid 
methods,  and  at  the  same  time  suggests  additive  correction  variants  that  lie  outside  the  scope  of  the  usual  mul¬ 
tigrid  formalism. 

For  example,  an  additive  correction  that  is  based  on  the  span  of  translates  of  "presmoothed"  iterates  does 
not  correspond  to  any  usual  "coarse  grid"  correction  in  the  multigrid  sense.  However,  it  can  be  shown  that  if 
such  corrections  follow  v  presmoothing  iterations,  and  if  the  smoothing  process  is  symmetrizable,  then  the  result 
is  at  least  as  good  as  that  which  would  have  been  obtained  by  doing  v  iterations  of  the  Chebyshev  semi-iterative 
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method  [A19]  with  optimum  parameters.  Furthermore,  a  model  problem  analysis  shows  that  these  corrections 
always  incorporate  the  ability  to  remove  completely  any  v  of  the  modes  that  remain  in  the  presmoothed  error. 
Details  of  these  results  will  appear  in  [31]. 

2.  The  Reduced  Basis  Method  for  Systems  of  Differential  Equations 

Error  Estimates  for  the  reduced  basis  method  solution  of  differential  and  differential-algebraic  equation 
systems  are  contained  in  the  PhJX  thesis  [2]  by  M.Y.  Lin-Lee,  under  the  direction  of  Professor  Thomas  Porsch- 
ing.  These  estimates  are  local  in  die  following  two  senses.  First,  they  apply  on  a  nontrivial,  but  possibly  very 
small  interval.  Second,  they  require  some  point  of  die  reduced  curve  to  lie  on  the  true  solution  curve.  The 
recent  research  reported  in  [18]  has  removed  the  interval  length  restriction  in  the  differential  equation  case  and 
extended  the  error  analysis  to  global  versions  of  the  methods  in  [2],  thus  effectively  eliminating  the  second  local 
aspect  of  that  work.  Furthermore,  this  work  also  incorporates  the  effects  of  the  errors  resulting  from  the  numeri¬ 
cal  integration  of  the  reduced  basis  ODE  systems. 

3.  Two-Fluid,  Two-Phase  Flow 

Additional  theoretical  results  on  the  nature  of  the  void  fraction  have  been  incorporated  into  [A6]  resulting 
in  the  revision  [3]. 

4.  Binary  Gas  Mixture  Flow  Through  Combustors 

In  an  attempt  to  reduce  the  development  cycle  costs  associated  with  design  of  gas  turbine  engine  combus¬ 
tors,  mathematical  combustor  models  are  being  employed  to  provide  information  about  performance  trends  and 
to  predict  velocity,  pressure  and  thermodynamic  property  profiles  in  simulated  practical  combustion  environ¬ 
ments.  It  has  been  demonstrated  that  the  dual  variable  method  can  be  applied  to  the  predictive  model  of  the 
fluid  dynamics  associated  with  an  axially  symmetric  centerbody  combustor  being  studied  at  WPAFB.  This  work 
appeared  in  [4],  Modifications  have  been  made  to  the  industry  standard  program  TEACH  to  allow  for  more 
efficient,  local  iteration.  The  same  algorithm  is  also  used  in  a  dual  variable  version  which  shows  even  more 
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savings. 


5.  Error  Estimation  and  Singularities 

The  central  theme  of  this  project  concerns  the  errors  arising  in  the  computational  solution  of  parameter- 
dependent  equations  of  the  form  F(zX)  =  0  where  F  is  a  given  nonlinear  mapping,  z  is  the  state  variable 
representing  die  solution,  aid  X  is  a  vector  of  parameters  that  characterize  intrinsic  properties  of  the  system  or 
extrinsic  quantities  influencing  its  behavior.  In  the  case  of  fluid  problems,  the  operator  may  be  generated,  for 
example,  by  the  time-independent  Navier  Stokes  equations  together  with  the  necessary  boundary  conditions. 

In  general,  the  solution  set  of  such  parametrized  equations  constitutes  a  differentiable  manifold  M  of 
dimension  equal  to  that  of  die  parameter  space  A.  While  this  fact  is,  of  course,  well  known,  we  appear  to  have 
been  the  first  who  have  been  using  this  fact  as  the  basis  for  a  consistent  study  of  the  numerical  problems  for 
these  equations.  Our  results  have  begun  to  show  clearly  the  value  and  power  of  this  geometric  approach. 

The  basic  procedures  for  the  computational  analysis  of  such  a  manifold  M  are  the  continuation  methods. 
When  M  has  dimension  p  >  1,  these  methods  require  a  restriction  to  some  path  on  M  and  then  produce  a 
sequence  of  points  along  that  path.  Obviously,  it  is  not  easy  to  develop  a  good  picture  of  a  multi-dimensional 
manifold  solely  from  information  along  such  paths.  This  led  us  recently  to  die  development  of  methods  for  the 
computation  of  simplicial  approximations  of  p-dimensional  subsets  of  M  (see  [25],  [26]).  Besides  these  methods 
for  computing  specific  sets  of  points  on  M,  another  important  class  of  numerical  procedures  concerns  the  detec¬ 
tion  and  computation  of  foldpdnts  on  M  with  respect  to  a  given  coordinate  space  T,  such  as  the  natural  parame¬ 
ter  space  A,  (see  [5],  [7],  [19],  [21]). 

For  numerical  considerations,  an  important  aspect  of  this  theory  concerns  the  estimation  of  the  various 
errors  arising  in  die  computations.  These  errors  fall  into  several  classes.  Fust  of  all,  there  is  the  important  ques¬ 
tion  of  the  discretization  errors  which  in  this  setting  turns  out  to  be  the  "distance"  between  the  manifolds  defined 
by  a  given  differential  equation  and  by  its  discretization,  respectively.  This  has  been  the  topic  of  a  series  of 
papers  by  J.  P.  Fink  and  W.  C.  Rheinboldt  with  partial  support  under  this  grant  See  the  earlier  articles  [A7], 
[A8]  and  then  [6]  and  [16].  In  particular,  in  the  last  named  paper  [16]  we  have  been  the  first  to  analyze  the  case 
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of  multi-dimensional  manifolds  which  is  of  increasing  importance  in  applications. 

An  essential  aspect  of  these  discretization  studies  concerns  the  question  as  to  the  exact  definition  of  the 
error  between  a  solution  manifold  and  its  approximation.  Obviously,  the  error  depends  on  which  points  are  to  be 
compared.  In  the  cited  papers  it  was  shown  that  this  correspondence  between  the  points  on  the  two  manifolds 
has  to  be  defined  by  appropriate  local  coordinate  systems.  In  other  words,  the  resulting  error  is  controlled  by  the 
choice  of  the  local  coordinate  system,  and,  since  the  error  measure  must  be  physically  meaningful,  not  all  local 
coordinate  systems  are  equally  desirable.  This  question  becomes  particularly  acute  in  the  vicinity  of  singular 
points  where  the  behavior  of  die  solutions  may  be  subject  to  change. 

This  connection  between  the  choice  of  the  local  coordinate  systems  and  the  singularity  behavior  of  a  point 
has  led  us  to  a  closer  study  of  admissible  coordinate  systems  at  foldpoints.  But  the  results  in  [5]  and  [7]  also 
suggest  that  the  open  questions  about  the  proper  choice  of  die  coordinate  system  for  physically  meaningful  error 
measures  requires  a  much  closer  look  at  the  characterization  of  foldpoints  on  manifolds  and  at  the  methods  for 
their  computation. 

Besides  die  discretization  errors  there  are  various  other  computational  errors  which  arise  in  the  analyse  of 
the  solution  manifolds  of  parametrized  equations.  In  particular,  both  in  the  continuation  methods  and  in  die  sim- 
pticial  approximation  methods,  mentioned  above,  we  encounter  die  following  sources  of  errors: 

(a)  Errors  associated  with  the  predictor  and  corrector 
(al)  Estimation  of  the  prediction  error 

(a2)  Control  of  die  termination  error 
(a3)  Control  of  the  effect  of  round-off 

(b)  Errors  associated  with  the  choice  of  local  coordinate 
system: 

(bl)  Estimation  of  the  domain  of  validity  of  this 
coordinate  system, 

(b2)  Control  of  the  range  of  applicability  of  the 
moving  frame  algorithm  in  the  simplicial 
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approximation  method. 

In  connection  with  foldpoints  computations,  various  further  error  questions  arise;  they  will  be  discussed  in  the 
next  subsection.  During  the  past  year  a  study  of  these  and  related  error  questions  has  begun  and  will  be  con¬ 
tinuing  during  die  next  year. 

6.  Detection  and  Computation  of  Singularities 

As  noted  in  the  previous  subsection,  the  solution  set  of  a  parametrized  equation  F(z,  X)  =  0  represents,  in 
general,  a  differentiable  manifold  in  the  combined  space  of  the  state  variable  and  the  parameter  vector.  This 
requires  a  regularity  assumption  which  is  not  very  restrictive  in  applications,  but  which  -  from  the  viewpoint  of 
singularity  theory  -  implies  the  use  of  a  suitable  "unfolding''.  In  most  practical  problems  a  host  of  very  natural 
unfoldings  suggest  themselves.  However,  the  particular  choice  of  unfolding  affects  the  location  and  type  of  the 
resulting  foldpoints  on  the  manifold  and  with  it  also  the  error  questions  raised  in  the  previous  subsection. 

Our  previous  work  in  this  area  concentrated  cm  a  study  of  the  differential  geometric  aspects  of  the  prob¬ 
lem.  In  the  already  mentioned  papers  [5],  [7]  we  used  the  tangent  map  to  develop  a  geometrically  instructive 
and  coordinate  free  treatment  of  foldpoints.  Then  in  [21]  is  was  shown  that  it  is  possible  to  reformulate  in  this 
setting  some  of  the  basic  results  of  bifurcation  theory  related  to  computational  aspects.  In  particular,  the  bifurca¬ 
tion  sets  appear  as  certain  cuts  of  the  solution  manifold.  This  in  turn  corresponds  to  the  consideration  of  particu¬ 
lar  augmented  equations  and  hence  opens  up  a  new  approach  to  the  study  of  effective  augmentations  for  die 
computation  of  singular  points.  Such  a  study  was  begun  in  [19], 

On  the  basis  of  the  results  in  [21]  R.  X.  Dai,  a  Ph.D.  student  at  the  University  of  Pittsburgh  under  the 
direction  of  Professor  Rheinboldt,  has  been  working  on  a  new  method  for  specific  class  of  foldpoint  problems. 
More  specifically,  a  new  minimal  augmentation  of  the  original  equations  have  been  devised  for  the  computation 
of  so-called  (p,l,l)  foldpoints.  For  the  resulting  systems  an  efficient  local  corrector  process  has  been  developed 
which  allows  for  the  continuation  along  a  path  of  such  foldpoints.  In  addition  the  approach  does  appear  to  open 
up  the  possibility  of  extending  the  simplicial  approximation  approach  mentioned  in  the  previous  subsection  to 
sub-manifolds  of  (p.1,1)  foldpoints.  This  is  currently  being  studied.  The  techniques  have  already  been  applied  to 
a  variety  of  test  problems  with  excellent  success.  A  comprehensive  report  will  be  forthcoming  within  a  few 
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months. 

As  noted  in  the  previous  subsection,  there  are  numerous  error  questions,  we  have  investigated,  arising  in 
the  different  processes  for  computing  foldpoints  on  a  solution  manifold.  We  mention  here  only  the  principal 
topics: 

(a)  Detea  a  foldpoint  in  a  neighborhood  of  a  given  set  x1 . x*  of  points. 

(b)  For  a  given  approximation  x  of  a  foldpoint  x*  of  the  manifold  estimate  the  distance  I  Ix-x*  1 1  between 
these  points. 

(c)  For  foldpoints  of  a  specific  type  develop  iterative  processes  which  are  locally  convergent  to  these  points 
from  any  sufficiently  close  approximation. 

(d)  If  x*  is  a  computed  foldpoint  of  the  discretized  equations  which  approximates  a  foldpoint  x  of  the  original 
differential  equations,  compute  an  estimate  of  the  distance  I  lx  -  x*  1 1. 

7.  Finite  Element  Formulation  of  the  Streanfuncdon-Vordcity  Equations 

The  Navier-Stokes  equations  can  be  written  in  primitive  variable  formulation,  in  terms  of  the  streamfunc- 
tion  as  a  fourth  order  problem  or  as  two  second  order  equations  in  the  streamfunction-vorticity  formulation.  In 
the  linear  case  the  fourth  order  problem  for  the  stream  function  is  the  well-known  biharmonic  equation. 
Although  the  primitive  variable  formulation  has  received  the  most  attention,  the  streamfunction-vorticity  formu¬ 
lation  is  also  of  considerable  interest  in  two  dimensional  domains.  That  is  partly  due  to  the  fact  that  only  two, 
as  opposed  to  three,  fields  are  to  be  computed;  but  it  is  mainly  due  to  the  fact  that  the  incompressibility  con¬ 
straint  is  avoided  through  the  introduction  of  the  streamfunction. 

Several  theoretical  and  practical  issues  arising  in  the  finite  element  approximation  of  the  streamfunction- 
vorticity  equations  have  been  studied.  Initially  error  estimates  for  the  linear  problem  were  investigated.  Since 
the  velocity  is  expressed  in  terms  of  the  derivatives  of  the  streamfunction,  it  is  of  practical  concern  to  ascertain 
if  these  derivatives  are  optimally  approximated  for  choices  of  elements.  Previous  analyses  concerning  this  prob¬ 
lem  were  improved  upon  and  the  optimality  of  the  error  verified  in  [A10], 
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Othcr  issues  arising  in  the  finite  dement  approximation  of  the  streamfunction  and  vorticity  include  compu¬ 
tations  in  multiply  connected  domains,  the  use  of  low  order  dements,  the  incorporation  of  a  variety  of  boundary 
conditions  into  the  weak  formulation,  estimates  for  the  errors  in  the  finite  approximations  for  the  nonlinear  prob¬ 
lem  and  the  recovery  of  the  primitive  variables.  A  preliminary  report  on  computations  in  multiply  connected 
domains  using  low  continuity  finite  element  spaces  was  presented  in  [11].  A  comprehensive  report  dealing  with 
all  of  the  theoretical  and  practical  issues  mentioned  above  as  well  as  presenting  numerical  examples  is  given  in 
[12]. 

8.  A  Finite  Element  Analysis  of  MHD  Flow 

The  equations  governing  the  steady  flow  of  incompressible  electrically  conducting  fluids  in  the  presence  of 
a  given  magnetic  field  can  be  expressed  as 

- Ait  +  i  (g  V  g)  +  -  <axV«)  -  Otxfixfi)  =  0 

n2  N 

+  div( m*£)  *  0 
div  m  =0 

where  m  is  the  velocity,  p  the  pressure,  $  the  electric  potential,  £  the  magnetic  field  and  N ,  n  given  dimension¬ 
less  parameters.  By  rewriting  certain  terms  using  vector  identities  and  using  appropriate  spaces,  one  can  obtain 
a  weak  formulation  for  this  problem  that  is  similar  to  one  for  the  Navier-Stokes  equations  written  in  terms  of 
primitive  variables  (see  [All]).  The  purpose  of  using  such  a  weak  formulation  is  to  take  advantage  of  the 
results  already  proved  in  [All].  Specifically,  the  weak  formulation  is  to  find  M  =  (g,<fr)  eff.pt  L,  such  that 

<*(M.k)  +  ai(g,iti)  +  f>(y,p)  =  (/ ,y)  for  ail  y  e  & 
h(M.X)-0  for  all 

where 

a(tf>  k)  =  -jp  I  Vu:Vy  +  J  {V^Ktfx£)}  {Vx|f-(yx£)} 
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J  (M'Vjtx  -«  V  vii] 

b(z,  x)  *  -/  X  X. 

andff  *a.‘  x/fj.L,1 2 3*  foel2:  /♦  =  OJ. 

a 

The  continuity  and  stability  conditions  necessary  to  guarantee  existence  and  uniqueness  of  the  solution  of 
the  weak  problem  have  been  proved.  In  addition,  an  error  estimate  for  the  finite  element  approximation  of  the 
weak  problem  has  been  obtained.  These  results,  as  well  as  a  discussion  of  an  iterative  method  for  solving  the 
discrete  problem  will  be  presented  in  [13]. 

9.  Dual  Variable  Transformations 

The  dual  variable  method,  originally  developed  in  the  context  of  finite  difference  discretizations  of  tran¬ 
sient  incompressible  Navier-Stokes  equations  [A9],  is  a  technique  to  considerably  reduce  the  size  of  the  linear 
system  to  be  solved  at  each  time  step.  The  method  involves 

(1)  the  determination  of  the  rank  of  the  discrete  divergence  operator,  A , 

(2)  the  determination  of  a  basis  for  the  null  space  of  A ,  and 

(3)  the  calculation  of  a  particular  solution  of  the  discrete  continuity  equation. 

In  [8]  a  finite  element  implementation  of  the  dual  variable  method  is  presented  using  quadrilateral  piece- 
wise  bilinear  velocity/constant  pressure  elements.  Algorithms  for  the  determination  of  a  basis  for  the  null  space 
of  the  discrete  divergence  operator  and  a  particular  solution  are  presented. 

In  [9]  a  finite  difference  discretization  of  the  Navier-Stokes  equations  describing  a  compressible  flow  prob¬ 
lem  is  viewed  as  a  system  defining  flows  on  an  associated  network.  This  observation  then  provides  a  means  of 
applying  the  dual  variable  method  to  economize  on  their  numerical  solution. 

The  nature  of  the  aerodynamics  in  and  around  such  structures  as  cavities  and  deflectors  or  spoilers  on  vari¬ 
ous  aircraft  configurations  was  investigated  using  the  dual  variable  method  [10]. 

A  summary  of  the  dual  variable  method  is  given  in  [14], 
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Iterative  methods  are  under  investigation  for  the  solution  of  the  linearized  finite  difference  discretizations 
involved  in  the  dual  variable  method.  The  generic  form  of  dual  variable  system  suggests  a  splitting  in  which  a 
Stieltjes  matrix  is  to  be  solved  at  each  step.  The  method  has  been  implemented  for  two  dimensional  domains 
and  convergence  properties  were  investigated  as  part  of  a  PhD.  thesis  by  George  Mesina  [35],  under  the  direc¬ 
tion  of  Professor  Charles  HalL 


10.  Fluid  Flaw  on  Curved  Domains 

A  finite  difference  scheme  was  derived  for  two-dimensional,  transient,  incompressible  Navier-Stokes  prob¬ 
lems  in  which  the  flow  domain  Q  is  a  bounded  simply-connected  region  for  which  there  exists  a  C2  invertible 
mapping  x  onto  the  unit  square: 

t :  Q  -*  S  ■  [0,l]x[0.1] 

The  transformed  Navier-Stokes  problem  is 

div  $  ■  4,>(  =  0 

$ j  n 

viJ  +  | j  j  viSj  =  ~Prjrj^i  +  |  j  |  (P/kv*^)ry  +  fi  1  =  1,2 

subject  to  initial  condition 

*(c.0)  =  a(*(c))  ,£  eS 


and  boundary  condition 


X.(L,  0  =  Q.£  edS  and  t  >  0  , 


where,  y(c,  t)  is  velocity,  p  is  pressure,  the  Jacobian  matrix 

J  * 

[Py]  *  »/  l/“1(/-1)r 

4=1/  l/‘,g  . 
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The  finite  difference  discretization  of  the  above  equations  is  proven  to  be  unconditionally  stable  and  con¬ 
vergent  They  also  reduce  to  the  well  known  Krzhivitsld  and  Ladyzhenskaya  scheme  [A12]  for  rectangular 
domains,  e.g.,  t  the  identity.  This  work  was  the  subject  of  a  Ph.D.  Dissertation  of  John  Ellison,  under  the 
direction  of  Professors  Charles  Hall  and  Thomas  Porsching  [A13]  and  appeared  in  [17]. 

11.  Differential  Equations  on  Manifolds 

As  noted  earlier,  equilibrium  problems  arise  naturally  in  continuum  and  fluid  mechanics  as  a  set  of  non¬ 
linear  equations  of  the  form 

£(*.&)  =  0  (1) 

where  i  is  from  the  state  space  Z  and  \  is  from  a  p  dimensional  parameter  space.  Then  the  equilibrium  mani¬ 
fold  is  the  set  of  solutions  to  (1)  in  D  cX  =  Zx A;  that  is 

M  -  {(i.WeS  I  £(*.&  =  (>}  (2) 

The  parameters  i  and  states  i  may  be  required  to  satisfy  a  differential  equation  of  the  form 

A(s)i  =  G(i)  (3) 

where  j  *  (i,  &).  We  then  interpret  (3)  as  a  differential  equation  on  the  manifold  (2),  DEM  for  short  Equa¬ 

tions  (1)  and  (3)  together, 

f  F  C*)  =  0 

=£(*) 

form  what  is  called  a  differential-algebraic  equation  (DAE). 

Two  applications  of  interest  to  the  investigators  in  which  DEM’s  occur  are: 

(i)  Incompressible  Fluid  Flow.  The  continuity  equation  is  an  algebraic  constraint  of  the  form  (1)  and  the 
time-dependent  Navier-Stokes  equations  are  the  differential  equation. 

(ii)  Punch  Stretching  of  Sheet  Metal.  The  principle  of  virtual  work  provides  a  force  equilibrium  equation 
which  defines  an  equilibrium  manifold  upon  which  one  seeks  solutions  to  differential  constitutive  laws  of 
the  form  (3). 
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In  [20],  a  new  numerical  method  was  presented  for  computer  simulations  of  punch  stretching  of  sheet 
metal.  Most  current  approaches  to  finite  element  modeling  of  large  deformation,  elastic-plastic  sheet  metal 
forming  use  a  rate  form  of  the  equilibrium  equations  and  then  must  correct  at  each  time  step  to  insure  that 
equilibrium  is  satisfied.  Such  methods  are  referred  to  as  incremental  methods.  The  new  method,  a  DEM 
approach  discretizes  the  more  fundamental  equilibrium  equations  in  non-rate  form  and  insures  equilibrium  of 
forces  at  each  time  step.  Formulating  the  problem  as  a  DEM  or  DAE  also  allowed  for  solution  of  the  discre¬ 
tized  system  using  off-the-shelf  software  such  as  LSODI.  Numerical  experimentation  indicated  that  the  DEM 
approach  was  computationally  much  more  efficient  than  the  incremental  approach. 

12.  Energy  Stability  of  Viscous  Incompressible  Flows 

The  problem  of  determining  sufficient  conditions  for  the  Sow  of  a  viscous  incompressible  fluid  to  be 
stable  under  arbitrary  disturbances  was  examined.  This  problem  is  of  importance  in  the  study  of  turbulence  and 
the  transition  which  occupies  a  region  of  space  and  is  subject  to  a  prescribed  velocity  distribution  on  the  boun¬ 
dary,  will  alter  radically  or  only  slightly  in  nature  if  it  is  perturbed  at  some  instant 

The  question  of  stability  can  be  addressed  by  either  standard  linear  perturbation  techniques  or  by  the 
energy  method;  the  latter  is  chosen  in  this  work.  Although  the  great  majority  of  stability  calculations  use  linear 
stability  analysis,  the  method  has  the  drawback  that  it  allows  only  perturbations  which  are  infinitesimal  in  mag¬ 
nitude.  This  rules  out  perturbations  of  finite  size  and  hence  cannot  give  accurate  information  in  many  cases. 
The  energy  method  allows  arbitrary  disturbances  but  its  shortcoming  is  that  the  disturbances  do  not  necessarily 
satisfy  the  Navier-Stokes  equations,  and  thus  the  stability  criterion  will  be  more  restrictive  than  in  the  actual 
physical  situation.  However,  the  energy  stsf  -lity  analysis  is  based  on  the  Navier-Stokes  equations  and  is  non¬ 
linear  in  nature  due  to  the  fact  that  no  linearizations  of  the  equations  are  done.  The  method  is  mathematically 
rigorous  and  does  give  insight  into  the  physical  situation. 

The  question  of  energy  stability  of  a  flow  can  be  formulated  as  a  linear  generalized  partial  differential 
eigenvalue  problem  even  though  the  analysis  is  based  on  the  nonlinear  Navier-Stokes  equations.  Essentially,  the 
procedure  is  to  obtain  an  equation  for  dK/dt  where  K  is  the  kinetic  energy  of  the  disturbance,  u,  and  then  to 
determine  conditions  which  guarantee  that  the  kinetic  energy  tends  to  zero  as  time  increases.  A  standard 
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eigen  value  problem  is  used,  where  stability  is  governed  by  the  dominant  eigenvalue.  Specifically,  we  have  that 
the  flow  is  stable  for  Reynolds  number  less  than  1 /X  where  X  is  the  dominant  eigenvalue  of  the  problem 

dg  -  grad  t  *  X  g  •  D 

div  g  *  0  ,  (1) 

g  =  0  on  the  boundary  . 

Here  D  is  the  deformation  tensor  of  the  unperturbed  flow. 

A  finite  element  method  is  used  to  approximate  the  dominant  eigenvalue  of  (1).  In  particular,  the  weak 
form  considered  is  to  find  nonzero  g  e  dj,  $  e  L2  and  X  e  R  such  that 

J  Vg:Vy  +  <i'v  g  -  X  J  g-D-g  for  all  y  e  dj 

(2) 

J  xdiv  g  =  0  for  all  %  e  L2 

where  L2  is  the  space  of  all  functions  which  are  square  integrable  and  tLl,  dJ  denote  the  usual  Sobolev  spaces. 
To  approximate  the  solution  to  the  weak  problem  (2),  finite  dimensional  subspaces  K*  c  dJ  and  Wh  c  L2  are 
chosen  which  depend  upon  a  parameter  0<A  <1  tending  to  zero.  The  approximate  problem  is  defined  analogous 
to  (2)  where  the  solution  is  sought  in  the  finite  dimensional  subspaces.  Once  bases  for  £*  and  Wh  are  chosen, 
the  approximate  problem  is  equivalent  to  an  algebraic  generalized  eigenvalue  problem.  An  estimate  for  the  error 
in  X  and  its  Galerkin  approximation  is  given  in  [A14], 

As  proposed,  a  code  was  developed  which  uses  a  mixed  finite  element  method  for  approximating  the  dom¬ 
inant  eigenvalue  of  (1).  The  program  was  used  to  determine  a  range  of  Reynolds  numbers  for  which  the  flow  is 
guaranteed  to  be  stable  for  the  examples  of  plane  shear  flow  and  Poiseuille  flow. 

The  first  example  is  the  simple  case  of  flow  in  a  channel  of  width  0<y<d  where  the  initial  velocity  is 
given  by  the  vector  (Ay  ,0),  the  deformation  tensor  is  given  by  D,j  =  0  for  i-j  and  Dij  =  ,5A  for  i  *  j, 
ij  =  1,2  and  the  Reynolds  number  is  kd2fh.  The  channel  is  assumed  infinite  in  length.  The  computed  Rey¬ 
nolds  number  for  a  channel  of  length  L  is  given  below. 

1/L  1  1/2  1/3  1/4  1/5 
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ReNo.  289.87  201.42  186.97S  182,11  180.01 


Using  the  above  calculations  the  extrapolated  value  at  1/1  =  0  is  177.4  which  is  in  good  agreement  with  the 
value  of  177  published  by  Orr. 


The  second  problem  of  determining  the  stability  of  Poise uille  flow  in  a  pipe  is  of  more  interest  physically 
and  has  been  studied  by  many  authors.  For  this  example  the  formulation  (1)  makes  sense  in  an  unbounded 
domain  when  the  solutions  are  single-valued  in  0  and  periodic  in  z.  To  this  end,  the  solution  is  assumed  of  the 
form  g(r,  0,  r)  =  l[(r,  a,  p)  where  P  is  an  integer.  This  form  is  substituted  into  (1)  and  the  follow¬ 

ing  system  results 


Lu  -  i 


2k 


-  -  -Ana 


Iv  +  i(2p/r2  u  -  P IrQ)  =  i) 


Lw  -i  a$  =  -Xru 


13,,. 

7  aT  *  ’ 


$. 


v  +  dw 


=  0 


V  J 

u,  v,  w  bounded  at  r  »  0 

«=v=w=0atr  =  l 


where  U  -  ( u,v,w ),  Lu  =  —  —  (ru)  - 

t  dr 


Bjfl  2 


u,  Lu  =  Lu  +  u/r2. 


Note  that  these  equations  form  a  complex  one-point  boundary  value  problem  with  a  singularity  at  the  ori¬ 


gin.  The  weak  formulation  must  incorporate  an  appropriate  boundary  condition  for  the  velocity  at  r  =  0.  The 
particular  weak  from  considered  is  the  following:  Find  u,  v,  w  e  Hu  $  e  H2,\eR  such  that  for  all  £  e  Hu 
Hi, 
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jur^r  dr  +  J  Cj-~  %dr  +  ij  vljdr  -  j^-^(r^jdr  =  \fwty2dr 

fv^rdr  +  J  Ci“§4r  -  iJ-22-uljdr  +  ijp^dr  =  0 

”  Jwr\,rdr  +  J  C2~§dr  +  ia|$ljrdr  *  kju%r2dr 

-J—(ri<)x^  -  i  J(lv  %dr  -  ijawxrdr  =  0 
where  Cx  ~  P2  +  1  +  cpr2  ,  C2=  Ct  -  1 

k 

for  appropriate  spaces  //t  aid  J/2- 

With  this  weak  formulation  the  condition  imposed  on  the  velocity  at  r  =  0  is  rur  =  rvr  =  nvr  =  0,  which 
is  a  natural  boundary  condition.  This  problem  was  discretized  using  piecewise  linear  elements.  The  results 
obtained  agree  with  those  of  Joseph  and  Carmi  [A15].  However,  their  numerical  calculations  were  unnecessarily 
complicated.  For  example,  different  techniques  for  the  various  cases  such  as  a  =  0  and  a  *  0  had  to  be 
employed  as  well  as  using  Frobenius  series  as  starting  values  near  the  origin.  Specifically,  the  value  of  81.5  was 
obtained  as  a  sure  limit  of  stability  and  corresponded  to  the  case  a  =  1,  P  *  1.  This  agreed  with  Joseph  and 
Carmi’s  result  and  confirmed  the  fact  that  the  value  of  R  *  180,  which  was  previously  believed  to  be  a  sure 
limit  of  stability,  is  incorrect.  This  work  is  discussed  in  a  paper  in  preparation  [23]. 

13.  The  Conjugate  Gradient  Method  on  Supercomputers. 

The  application  of  die  preconditioned  conjugate  gradient  (PCCG)  method  to  the  solution  of  the  large, 
sparse  linear  systems  resulting  from  die  discretization  of  partial  differential  equations  on  regular  and  irregular 
domains  was  investigated.  The  primary  goal  was  the  efficient  implementation  of  the  PCCG  method  on  vector 
supercomputers.  In  [24],  this  goal  was  met  by  the  introduction  of 

(i)  A  data  structure  which  is  suitable  for  manipulation  on  vector  machines, 

(ii)  Preconditioned  matrices  which  are  effective  for  general  sparse  matrices,  and 

(iii)  Numbering  schemes  for  both  regular  and  irregular  grids,  which  are  amenable  to  adequate  vectorizadons. 

14.  The  Numerical  Solution  of  Transport  Problems 

The  question  of  finding  schemes  for  approximating  the  solution  of  transport  problems  (e  =  0)  and 
diffusive  transport  problems  (0  <  e  «  0(1))  is  a  central  problem  in  computational  fluid  dynamics  and  combusdon 
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theory.  Frequently,  efficient  schemes  are  required  which  have  high  accuracy,  preserve  positivity,  minimize 
dispersive  effects,  minimize  grid  orientation  effects,  etc. 

An  additive  correction  algorithm  proposed  in  [A17]  was  considered  for  the  1-D  equilibrium  problem: 


f 


l-e  u"  +  f  (x)u  +  g(x)u  =  gQc)  , 0 <x  <  1 
u(0)  =  a  ,  u(l)  =  P  . 


(1) 


An  exact  analysis  of  the  deferred  correction  algorithm  was  performed  in  [27],  accounting  for  layer  effects,  via 
the  theory  of  discrete  barrier  functions.  This  analysis  established  that  in  outer  layer  regions  (regions  of  fixed 
distance  5  >  0  from  layers)  the  scheme  converges  uniformly  in  e  with  accuracy  given  by  the  high  accuracy 
"defect"  operator.  The  influence  of  discrete  layer  terms  was  clarified  in  [27]  also.  A  mesh  refinement  process  is 
naturally  suggested  by  this  analysis,  which  also  gives  a  bound  upon  the  region  of  assured  Offi2)  accuracy. 

Next  the  case  of  2-D  equilibrium  transport  (e  =  0)  and  diffusive  transport  (e  >  0)  problems  was  considered 
u  =  u(x,  y )  satisfies 


t 


j-e  An  +  Vi  h,  +  v2«,  +  gu  =  q  ,  in  Q  c  R2  , 
ji  =  a  on  XI ,  0^e<  0(1)  . 


(2) 


We  have  focused  our  intention  on  methods  which  strictly  preserve  positivity.  [A  2-D  filtering  technique  which 
approximately  preserves  positivity  in  die  deferred  correction  methods  previously  discussed  was  introduced  in 
[A18]].  The  simplest  way  to  satisfy  this  requirement  is  with  difference  stencils  which  are  of  positive  type. 

In  [28]  the  spurious  diffusive  effects  of  positive  type  schemes  were  investigated  for  (2).  Specifically,  the 
presence  of  anisotropic  crosswind  diffusion  was  studied.  It  was  shown  in  [28]  that,  except  when  the  transport 
direction  is  precisely  aligned  with  the  mesh,  there  does  not  exist  a  positive  type  scheme  with  zero  spurious 
crosswind  diffusion.  The  precise  amount  of  this  diffusive  term  is  given  in  [28].  In  [32],  this  analysis  was  taken 
one  step  further.  The  principal  axes  of  diffusion  for  difference  schemes  woe  defined  By  comparing  these  to 
the  physical  flow,  a  general  methodology  for  comparing  difference  schemes  and  parameter  selection  in  a  given 
scheme  is  obtained 


In  parallel,  a  much  more  ambitious  approach  to  2-D  transport  problems  has  been  investigated.  Monotone 
(i.e.  inverse  positive)  type  schemes  for  transport,  transport-diffusion  and  transport-diffusion-reaction  equations 
have  been  studied.  Monotone  schemes  can,  in  principle,  overcome  many  of  the  inherent  barriers  in  positive  type 
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methods.  However,  the  inherent  global  and  nonlinear  character  in  generating  monotone  stencils  presents  a 
severe  challenge. 

Recently,  there  have  been  three  breakthroughs  at  ICMA  on  this  question.  The  first  is  that  a  methodology 
has  been  developed  for  reducing  to  sequential  linear  systems  these  global  nonlinear  consistency  conditions. 
Thus,  a  technique  has  been  developed  for  generating  a  class  of  novel  monotone  type  difference  schemes  for  gen¬ 
eral  boundary  value  problems.  This  class  contains  a  number  of  free  parameters  which  can  be  used  to  enforce 
other  desirable  features  e.g.,  high  accuracy,  minimal  dispersive  or  diffusive  properties,  etc. 

Secondly,  a  "collage"  type  theorem  for  monotone  and  positive  type  matrices  has  been  proven.  General 
monotone  (inverse  positive)  matrices  lack  the  algebraic  structure  of  matrices  of  positive  type.  Thus,  it  is  neces¬ 
sary  to  give  a  systematic  procedure  for  combining  monotone  type  stencils  with  other  stencils  to  preserve  stability 
in  complex  convection-diffusion-reaction  problems.  This  collage  theorem  is  a  major  step  to  developing  a  gen¬ 
eral  theory  which  will  specify  the  corcect  procedures. 

These  breakthroughs  are  a  very  general  approach/methodology  for  deriving  schemes  with  assured  stability 
and  assured  preservation  of  positivity.  To  complete  the  picture,  we  are  currently  working  on  generating  large 
classes  of  novel  schemes  for  2-D  transport  problems  and  analysis  of  the  error  the  schemes  generated.  These 
would  then  be  applied  to  many  specific  problems  including  combustion  and  fluid  flow. 

15.  Problems  with  Memory  in  Mathematical  Physics 

Progress  has  been  made  in  a  tangentially  related  area  of  applied  mathematics:  physical  processes  incor¬ 
porating  delay  on  memory  effects.  For  problems  with  (possibly  infinite)  delay: 

► 

x(t)+  J  g(t,s,x(t),x(t-s))d\xis)  =  /(/),/  e  R*  , 

x(s)  =  #s)  ,  je/r, 

k 

a  global  existence  theory  has  been  developed  in  [29]  where  the  asymptotic  behavior  of  solutions  is  also  con¬ 
sidered.  In  a  companion  paper  [30],  general  conditions  upon  g  and  p  are  given  which  ensure  exponential 
asymptotic  stability  of  the  initial  value  problem.  Thus,  [29],  [30]  give  a  fairly  complete  picture  of  the  behavior 
of  solutions  to  (1)  under  a  very  mild  nonresonance  type  condition  upon  g . 
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16.  3-D  Fluid  Flow  Calculations  on  the  CRAY  XMP-48 

There  is  a  rather  large  (6400  flow  cell)  version  of  the  2-D  fluids  code  ALGAE  [10]  which  typically  takes 
overnight  on  a  VAX  8650  to  get  one  time- step  finished.  We  needed  such  a  large  version  in  a  fluids  simulation 
of  aircraft  cavity  flows,  and  found  that  we  could  get  300  time  steps  (often  needed  to  get  a  steady  state)  finished 
on  a  CRAY  XMP-48  within  a  24  hour  period.  However,  we  noticed  that  this  code  was  extremely  slow  and  that 
75-90%  of  the  CPU  time  on  the  CRAY  was  being  used  by  the  solver.  One  obvious  conclusion  was  that  any 
practical  3-D  code  would  have  to  be  developed  in  a  supercomputing  environment  using  very  efficient  vectorized 
solvers. 

Because  of  a  short  term  need  to  make  die  CRAY  XMP-48  version  of  ALGAE  more  efficient,  and  to  give 
us  some  experience  with  high  efficiency  solvers  of  linear  systems  on  the  CRAY,  we  attempted  to  unravel  the 
bookeeping  in  ALGAE  to  produce  the  final  equations  in  matrix-vector  form  without  turning  on  ALGAE’s 
inefficient  (on  the  CRAY)  frontal  solver.  This  latter  feature,  combined  with  ALGAE’s  link-list  bookeeping,  pro¬ 
vides  its  excellent  capability  for  efficient  solution  of  a  complicated  geometry  on  computers  with  limited  storage. 
However,  on  a  large-memary  machine  devoted  to  vector  efficiency  and  with  parallel  capability,  the  frontal  solu¬ 
tion  can  be  hard  to  vectorize,  and  link-list  bookeeping  requires  re-organization  in  vector  form  to  get  reasonable 
efficiency  on  a  supercomputer  such  as  the  XMP-48.  We  had  hoped  to  modify  ALGAE  to  simply  generate  the 
matrices  without  turning  on  the  frontal  solver  so  we  could  call  one  of  the  high  efficiency  solver  packages  on  the 
CRAY. 

This  turned  out  to  be  a  coding  nightmare,  so  it  was  decided  to  move  on  to  a  new  coding  for  the  3-D  ver¬ 
sion.  We  have  extended  the  dual  variable  method  [A9]  to  incompressible  flow  in  three  dimensions  and  a  gradu¬ 
ate  student,  Ms.  Ye,  under  the  direction  of  Professor  Charles  Hall  has  developed  an  algorithm  for  the  construc¬ 
tion  of  "countries”  at  faces  of  the  associated  network.  A  node-oriented  system  has  been  planned  on  domains  that 
are  imbedded  in  rectangular  parallelepipeds. 
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